//This code is an example for calculting the hardcore boson dynamics from finite T initial states
//Output file is the one-body correlation matrix <b^\dagger_i b_j>
//Ref: DOI: 10.1103/PhysRevA.95.033617

#include <complex>
#define MKL_Complex16 std::complex<double>
#include <iostream>
#include <fstream>
#include <mkl.h> 
std::complex <double> i_imag(0, 1);
using namespace std;
double Z, Z_coeff;

double *Matrix(double V, double t, int N);
double H_0(int braNumber, int ketNumber, double V, double t, int N);
double Find_mu(double* Energylevels, int Hdim, int N_Expected, double T, double mu_1, double mu_2, double delta_mu);
std::complex <double>* GE_dt(int N, double* U, double* E, double dt);
std::complex <double>* Umulti(int N, std::complex <double>* U1, std::complex <double>* U2);
std::complex<double>* StateHCB(std::complex<double>* U, double* Energylevels, int Hdim, double T, double mu);
std::complex <double> Gmn_GE_temp(int m, int n, int N, std::complex<double>* mat_temp);
std::complex <double> determinant_complex_Z(std::complex <double>* mat, int N);

int main(int argc, char * argv[])
{
	int N = 50;				    //Chain size
	int Nk = 5;				    //Particle number
	double time = 5;            //Evolution time
	double V_init = 1e-2;		//Init Harmonic trap
	double V_final = 5e-2;		//Final Harmonic trap     
	double t_hopping = -1;		
	double temperature = 0.1;   //Temperature

	double *Hmatrix_init = NULL;
	int info_init;
	double *w_init = NULL;
	Hmatrix_init = Matrix(V_init, t_hopping, N);
	w_init = new double[N];
	info_init = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'U', N, Hmatrix_init, N, w_init);
	if (info_init != 0)
	{
		cerr << info_init << "\n"; exit(0);
	}
	double *Hmatrix_fin = NULL;
	double *w_fin = NULL;
	int info_fin;
	Hmatrix_fin = Matrix(V_final, t_hopping, N);
	w_fin = new double[N];
	info_fin = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'U', N, Hmatrix_fin, N, w_fin);
	if (info_fin != 0)
	{
		cerr << info_fin << "\n"; exit(0);
	}
	double mu;
	mu = Find_mu(w_init, N, Nk, temperature, -1000, 1000, 1e-10);
	Z = 1; Z_coeff = 0;
	for (int i = 0; i < N; i++)
	{
		Z = Z * (1 + exp(-(w_init[i] - mu) / temperature));
		while (abs(Z) > 10)
		{
			Z = Z / 10;
			Z_coeff++;
		}
		while (abs(Z) < 1)
		{
			Z = Z * 10;
			Z_coeff--;
		}
	}
	std::complex<double> *Utemp0 = NULL;
	std::complex<double> *Utemp1 = new std::complex<double>[N*N];
	std::complex<double> *Utemp2 = NULL;
	std::complex<double>* state = NULL;
	Utemp0 = GE_dt(N, Hmatrix_fin, w_fin, -time);
	for (int ii = 0; ii < N; ii++)
	{
		for (int jj = 0; jj < N; jj++)
		{
			Utemp1[ii*N + jj] = Hmatrix_init[jj*N + ii];
		}
	}
	Utemp2 = Umulti(N, Utemp1, Utemp0);
	for (int ii = 0; ii < N; ii++)
	{
		for (int jj = 0; jj < N; jj++)
		{
			Utemp0[ii*N + jj] = conj(Utemp2[jj*N + ii]);
		}
	}
	delete[] Utemp2;
	state = StateHCB(Utemp0, w_init, N, temperature, mu);
	std::complex <double>* corrmat = NULL;
	std::complex <double> resultall;
	corrmat = new std::complex <double>[N*N];
	for (int i0 = 0; i0 < N; i0++)
	{
		for (int j0 = i0; j0 < N; j0++)
		{
			resultall = Gmn_GE_temp(i0, j0, N, state);
			corrmat[i0*N + j0] = resultall;
			if (i0 != j0)
			{
				corrmat[j0*N + i0] = std::conj(resultall);
			}
		}
	}
	std::ofstream outputdata;
	outputdata.open("Corrmat.dat", ios::out | ios::binary);
	outputdata.write((char *)corrmat, sizeof(std::complex<double>)*N*N);
	outputdata.close();
	delete[] Hmatrix_init;
	delete[] w_init;
	delete[] Hmatrix_fin;
	delete[] w_fin;
	delete[] corrmat;
	delete[] Utemp0;
	delete[] Utemp1;
	delete[] state;
}


double *Matrix(double V, double t, int N)
{
	int bra, ket;
    double *matrix=new double [N*N];
	for (bra=0;bra< N; bra++)
    {
       for (ket=0; ket< N; ket++)
	   {
		   matrix[bra*N+ket] = H_0(bra, ket, V, t, N); 
       }
    }
    return matrix;
}

double H_0(int braNumber, int ketNumber, double V, double t, int N)
{
	double energy = 0;
	int temp, i, n;
	if ((abs(braNumber) == abs(ketNumber)))
	{
		energy = energy + V * (braNumber - N / 2)*(braNumber - N / 2);
	}
	if ((braNumber==ketNumber+1)||(braNumber==ketNumber-1))
	{
		energy = energy + t;
	}
	return energy;
}



double Find_mu(double* Energylevels, int Hdim, int N_Expected, double T, double mu_1, double mu_2, double delta_mu)
{
	int i;
	double N_1, N_2, N_test;
	double mu;
	double delta; 
	N_1 = 0;
	N_2 = 0;
	N_test = 0;
	for (i = 0; i < Hdim; i++)
	{
		N_1 += 1 / (1 + exp((Energylevels[i] - mu_1) / T));
		N_2 += 1 / (1 + exp((Energylevels[i] - mu_2) / T));
	}
	mu = (mu_1 + mu_2) / 2;
	if ((N_1 < N_Expected) && (N_2 > N_Expected))
	{
		while (mu_2 - mu_1 > delta_mu)
		{
			mu = (mu_1 + mu_2) / 2;
			for (i = 0; i < Hdim; i++)
			{
				N_test += 1 / (1 + exp((Energylevels[i] - mu) / T));
			}
			if (N_test > N_Expected)
			{
				mu_2 = mu;
			}
			else if (N_test < N_Expected)
			{
				mu_1 = mu;
			}
			else
			{
				return mu;
			}
			N_test = 0;
		}
		return mu;
	}
	else if ((N_1 > N_Expected) && (N_2 < N_Expected))
	{
		while (mu_2 - mu_1 > delta_mu)
		{
			mu = (mu_1 + mu_2) / 2;
			for (i = 0; i < Hdim; i++)
			{
				N_test += 1 / (1 + exp((Energylevels[i] - mu) / T));
			}
			if (N_test > N_Expected)
			{
				mu_1 = mu;
			}
			else if (N_test < N_Expected)
			{
				mu_2 = mu;
			}
			else
			{
				return mu;
			}
			N_test = 0;
		}
		return mu;
	}
	else
	{
		if (N_1 == N_Expected)
		{
			return mu_1;
		}
		else if (N_2 == N_Expected)
		{
			return mu_2;
		}
		else
		{
			cout << "Initial mu guass error \n";
			exit(0);
		}
	}
}

std::complex <double>* GE_dt(int N, double* U, double* E, double dt)
{
	std::complex <double>* result_temp = new std::complex <double>[N*N];
	std::complex <double>* result = new std::complex <double>[N*N];
	std::complex <double>* Ucomplex = new std::complex <double>[N*N];
	int i, j;
	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			result_temp[i*N + j] = exp(-i_imag * E[i] * dt)*U[j*N + i];
			Ucomplex[i*N + j] = U[i*N + j];
		}
	}
	std::complex <double> alpha(1, 0);
	std::complex <double> beta(0, 0);
	cblas_zgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,N,N,N,&alpha,Ucomplex,N,result_temp,N,&beta,result,N);
	delete[] result_temp;
	delete[] Ucomplex;
	return result;
}

std::complex <double>* Umulti(int N, std::complex <double>* U1, std::complex <double>* U2)
{
	std::complex <double>* result = new std::complex <double>[N*N];
	std::complex <double> alpha(1, 0);
	std::complex <double> beta(0, 0);
	cblas_zgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,N,N,N,&alpha,U1,N,U2,N,&beta,result,N);
	return result;
}


std::complex<double>* StateHCB(std::complex<double>* U, double* Energylevels, int Hdim, double T, double mu)
{
	long dim;
	dim = Hdim * Hdim;
	std::complex<double> *corr_mat = new std::complex<double>[dim];
	int m, n;
	std::complex<double> *temp_mat = new std::complex<double>[dim];
	std::complex<double> *U_dagger = new std::complex<double>[dim];
	double Inverse_temp;
	for (m = 0; m < Hdim; m++)
	{
		Inverse_temp = exp(-(Energylevels[m] - mu) / T);
		for (n = 0; n < Hdim; n++)
		{
			U_dagger[m*Hdim + n] = conj(U[n*Hdim + m]) * Inverse_temp;
		}
	}
	//matrix mulplication
	std::complex <double> alpha(1, 0);
	std::complex <double> beta(0, 0);
	cblas_zgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,Hdim,Hdim,Hdim,&alpha,U,Hdim,U_dagger,Hdim,&beta,temp_mat,Hdim);
	delete[] U_dagger;
	return temp_mat;
}

std::complex <double> Gmn_GE_temp(int m, int n, int N, std::complex<double>* mat_temp)
{
	int i, j;
	std::complex <double>* mat0 = NULL;
	std::complex <double>* mat1 = NULL;
	std::complex <double> result;
	std::complex <double> det0, det1;
	if (m == n)
	{
		mat0 = new std::complex <double>[N*N];
		mat1 = new std::complex <double>[N*N];
		for (i = 0; i < N; i++)
		{
			for (j = 0; j < N; j++)
			{
				if (i < n && j >= m)
				{
					mat0[i*N + j] = -mat_temp[i*N + j];
				}
				else if (i >= n && j < m)
				{
					mat0[i*N + j] = -mat_temp[i*N + j];
				}
				else
				{
					mat0[i*N + j] = mat_temp[i*N + j];
				}
			}
		}
		for (i = 0; i < N; i++)
		{
			for (j = 0; j < N; j++)
			{
				if (i == m)
				{
					mat1[i*N + j] = -mat0[i*N + j];
				}
				else
				{
					mat1[i*N + j] = mat0[i*N + j];
				}
			}
		}
		for (i = 0; i < N; i++)
		{
			mat0[i*N + i] += 1;
			mat1[i*N + i] += 1;
		}
		det0 = determinant_complex_Z(mat0, N);
		det1 = determinant_complex_Z(mat1, N);
		result = (-det1 + det0) / (2);
		delete[] mat0;
		delete[] mat1;
		return result;
	}
	else
	{
		mat0 = new std::complex <double>[N*N];
		mat1 = new std::complex <double>[N*N];
		for (i = 0; i < N; i++)
		{
			for (j = 0; j < N; j++)
			{
				if (i < n && j >= m)
				{
					mat0[i*N + j] = -mat_temp[i*N + j];
				}
				else if (i >= n && j < m)
				{
					mat0[i*N + j] = -mat_temp[i*N + j];
				}
				else
				{
					mat0[i*N + j] = mat_temp[i*N + j];
				}
			}
		}
		for (i = 0; i < N; i++)
		{
			for (j = 0; j < N; j++)
			{
				mat1[i*N + j] = mat0[i*N + j];
				if (i == m)
				{
					mat1[i*N + j] += mat0[n*N + j];
				}
			}
		}
		for (i = 0; i < N; i++)
		{
			mat0[i*N + i] += 1;
			mat1[i*N + i] += 1;
		}
		det0 = determinant_complex_Z(mat0, N);
		det1 = determinant_complex_Z(mat1, N);
		result = (det1 - det0);
		delete[] mat0;
		delete[] mat1;
		return result;
	}
}

std::complex <double> determinant_complex_Z(std::complex <double>* mat, int N)
{
	int i, info;
	int *ipiv = new int[N];
	std::complex <double> det1, result;
	det1 = 1 + 0 * i_imag;
	double coeff = 0;
	info = LAPACKE_zgetrf(LAPACK_ROW_MAJOR, N, N, mat, N, ipiv);
	for (i = 0; i < N; i++)
	{
		if (ipiv[i] != (i + 1))
		{
			det1 = -det1;
		}
		det1 = det1 * mat[i*N + i];
		while (abs(det1) > 10)
		{
			det1 = det1 / 10;
			coeff++;
		}
		while (abs(det1) < 1)
		{
			if (det1 == 0) { result = 0; return result;}
			det1 = det1 * 10;
			coeff--;
		}
	}
	det1 = det1 / Z;
	coeff = coeff - Z_coeff;
	result = det1 * pow(10, coeff);
	delete[] ipiv;
	return result;
}